#include "stdafx.h" 
#include "sequence.h" 
#include "timer.h" 
#include "strutils.h"


Sequence::~Sequence ()//'~' is a destructor, which we use to clean up our memory. Every Malloc must have a free here
{
	free (sequenceData); 
	free (hashTable); //whenever we create pointers to memory you need to set them zero in the construtor and free them in the destructor 
	free (locationsInLibrary);
}

Sequence::Sequence ()
{
	//this is the constructor, it gets called automatically when we "instantiate an instance of sequence"
	//initialize member data
	init();
}

unsigned int Sequence::ReadFromFile(ifstream & seqin)
{
		SequenceNumber = 0; 
	memset (ascensionNumber, 0, sizeof (ascensionNumber)); //fill every byte with zeroes up to the size of the ascentionNumber   
	memset (Title, 0, sizeof (Title)); //fill every byte with zeroes up to the size of the Title  
	sequenceData = 0;
	const int bufferSize = 1000; //"buffer" is our memory scratch pad, needs a placeholder to put the data from the file 
	char buffer[bufferSize] ={0}; //fill with all zeroes 
	sequenceData = (char*) malloc (2000000); //malloc returns a generic address (void*) so we need to tell its a char* 
	memset (sequenceData, 0, 2000000); 
	//read the header
	seqin.getline(header, sizeof (header), '\n'); //read until I hit carriage return 

	//read until i hit next header 
	while(!seqin.eof())//keep going until I reach end of file 
	{

		seqin.getline(buffer, bufferSize, '\n'); //read file, up to bufferSize bytes until you hit a carriage return

		if (buffer[0]== '>')
		{

			break;
		}
		else
		{
			strcat(sequenceData,buffer);
			//cout << buffer << "\n"; 
		}
	}
	sequenceLength = strlen (sequenceData); 
	return sequenceLength;
}

void Sequence::init(unsigned int librarySize,unsigned int numberOfNucleotidesInHash,  unsigned int randomSeed)
{
	SequenceNumber = 0; 
	header[0]=0;
	expectedHashCollisions = -1;
	memset (ascensionNumber, 0, sizeof (ascensionNumber)); //fill every byte with zeroes up to the size of the ascentionNumber   
	memset (Title, 0, sizeof (Title)); //fill every byte with zeroes up to the size of the Title  
	sequenceData = 0;
	sequenceLength = 0;
	NucleotiesInHash = numberOfNucleotidesInHash; 
	hashTable = 0; 
	locationsInLibrary = 0;
	hashCollisions=-1; 
	int expectedHashCollisions=-1;
	timeToCreateHash=0.0;
	sequenceLength = librarySize;
	sequenceData = (char*) malloc (librarySize+1); //malloc returns a generic address (void*) so we need to tell its a char* 
	memset (sequenceData, 0, librarySize+1); 
	if(randomSeed!=0)
	{
		//fill with random data based on seed
		PopulateRandomSequence(randomSeed,librarySize, sequenceData);
	}

}


Sequence::Sequence (unsigned int librarySize, unsigned int numberOfNucleotidesInHash, unsigned int randomSeed )
{
	init( librarySize, numberOfNucleotidesInHash,randomSeed );
	

}
	

Sequence::Sequence (const char* filename)
{
	//this is the constructor, it gets called automatically when we "instantiate an instance of sequence"
	//initialize member data
	init();
	ifstream in (filename); 
	ReadFromFile (in); 
}



bool Sequence::allocateHashTable()
{
	if(hashTable != 0)
	{
		free (hashTable);
		hashTable=0;
	}
	hashTable = (int*)malloc (getHashTableSize() * sizeof(int)); //need enought room for hash table length intergers, getting memory (malloc) 
	memset (hashTable, 0, getHashTableSize() * sizeof(int));
	return (true); 
}

int GetHashValue(char * pLib, int offset, int nucleotidesPerHash)
{
		int HashVal =0;
		for (int j = 0; j < nucleotidesPerHash; j++)
		{
			char Nucleotide = pLib [offset+j];
			int NucValue = 0; 
			switch (Nucleotide)
			{
			case 'A': 
				{
					NucValue = 0;
				}
				break;
			case 'G': 
				{
					NucValue = 1;
				}
				break;
			case 'T': 
				{
					NucValue = 2;
				}
				break;
			case 'C': 
				{
					NucValue = 3;
				}
				break;
			default:
				break;
			}

			HashVal = HashVal << 2 ; //The '<<' operator will shift bits to the left when the char on the left is a number and not a 'cout' 
			HashVal += NucValue;
		}//end of Jloop

		return HashVal ;
}

int GetHashValue(unsigned int * pLib, int offsetBytes, int nucleotidesPerHash)
{
	//This hash function assumes that the characters are encoded four each into an int array
	//  Offest is still in chars
	//  int1   int2  int3
	//[ ACGT] [CGTA] [GTAC]

	//First, we need to figure out how walk bytes in the stream starting at offset
	// Divide by 4 will get the correct integer -- so for example, if we start at byte 2, 2/4 = 0 in integer division
	// we use the modulo (remainder to figure out how much to shift
	

	
		int HashVal =0;
		for (int j = 0; j < nucleotidesPerHash; j++)
		{
		// Divide by 4 will get the correct integer -- so for example, if we start at byte 2, 2/4 = 0 in integer division
		unsigned int cval =  pLib [j/4];
		// we use the modulo (remainder to figure out how much to shift
		int remainder = j % 4;

		//Integers are stored in backwards -- least significant byte first  (ACGT in memory corresponds to TGCA as a 32 bit integer)
		// if the remainder is 0, we are at 0, 4, 8, etc    ACGT so we want A
		//			which encodes least significant byte  so we need to shift to right 3*8 bits
		// if the remainder is 1, we are at 1, 5, 9, etc    ACGT so we want C,
		//			so we need to shift to right 1*8 bits
		// if the remainder is 2, we are at 2, 6, 10, etc   ACGT so we want G, so we need to shift to right 2*8 bits
		// if the remainder is 3, we are at 3, 7, 11, etc   ACGT so we want T, so we need to shift to right 3*8 bits
		cval = cval >>  remainder *8;
	
		//now we need to mask off the upper three bytes by ANDing with 0x000000FF;
		cval =cval & 0x000000FF;

	
			int NucValue = 0; 
			switch (cval)
			{
			case 'A': 
				{
					NucValue = 0;
				}
				break;
			case 'G': 
				{
					NucValue = 1;
				}
				break;
			case 'T': 
				{
					NucValue = 2;
				}
				break;
			case 'C': 
				{
					NucValue = 3;
				}
				break;
			default:
				break;
			}

			HashVal = HashVal << 2 ; //The '<<' operator will shift bits to the left when the char on the left is a number and not a 'cout' 
			HashVal += NucValue;
		}//end of Jloop

		return HashVal ;
}


int Sequence::LocateSequenceInLibrary(Sequence & library)
{
	free(locationsInLibrary);
	 locationsInLibrary  = (int*) malloc (sizeof(int) * sequenceLength );
	 for (int seqLoc = 0; seqLoc < (sequenceLength - NucleotiesInHash); seqLoc++)
	 {
		int HashVal =GetHashValue(sequenceData,seqLoc,NucleotiesInHash);
		locationsInLibrary[seqLoc]=library.hashTable[HashVal]-seqLoc;
	 }
	 return 0;
}

int Sequence::ResolveAlignment(Sequence & library)
{
	if(locationsInLibrary == NULL)
		LocateSequenceInLibrary(library);

	std::sort(&locationsInLibrary[0], &locationsInLibrary[sequenceLength-NucleotiesInHash-1]);

	int commonValue = 0; 
	int commonValueCount = 0;
	int maxCommonValueCount =0;
	int maxCommonValue =0;
	int zeroCounter = 0; 

	for (int i = 0; i < sequenceLength - NucleotiesInHash; i++)
	{
		int value = locationsInLibrary[i] ;
		if (value == 0)
		{
			zeroCounter++; 
			continue; 
		}
		if (commonValue == value)
		{
			commonValueCount++;
			if (commonValueCount > maxCommonValueCount)
			{
				maxCommonValueCount = commonValueCount;
				maxCommonValue = value; 
			}
		} 
		else
		{
			if (commonValueCount > 2)
			{
				cout << "Found " << commonValue << " a total of " << commonValueCount <<  " times\n";
			}
			commonValue = value; 
			commonValueCount = 0; 
		}
	}
	
	//summarize results:
	int numberOfPotentialHits = sequenceLength - NucleotiesInHash;
	int numberOfActualHits = numberOfPotentialHits - zeroCounter;
	double HashCollisionRate = (double) library.expectedHashCollisions / (double) library.sequenceLength;
	double libraryFillRate =  ((double) library.sequenceLength - (double) library.expectedHashCollisions )/ (double) library.getHashTableSize();
	int numberOfExpectedCommonMode =numberOfPotentialHits - (numberOfPotentialHits *  HashCollisionRate );
	int numberOfExpectedRandomHits = numberOfPotentialHits * libraryFillRate;
	double likelihoodOfFindingHitsonAlltries=pow(libraryFillRate,numberOfPotentialHits);
	double likelihoodOfFindingNhitsCommonMode = pow(NucleotiesInHash + maxCommonValueCount,4.0);

	
	cout << "------Results for Probe Search " << this->Title << " in Library "  << library.Title << "\n";
	cout << "\t\tLikelihood  of Hitting on All Hashees is  1 in " << (1.0/likelihoodOfFindingHitsonAlltries) << ")\n";
	cout << "\t\tHIT ON " << numberOfActualHits << " out of " <<   numberOfPotentialHits << " possible hash table hits\n";
	cout <<  "\t\tExpected common mode count based on collision rate is " << numberOfExpectedCommonMode << "\n";
	cout << "\t\tFound a common mode of " << maxCommonValue << " a total of "  << maxCommonValueCount << " times \n";
	cout << "\t\tLikelihood  of finding this many common mode hits is 1 in " << likelihoodOfFindingNhitsCommonMode << "\n";
	cout << "\t\tCONCLUSION: ";
	if(sequenceLength / likelihoodOfFindingNhitsCommonMode < .0001)
	{
		cout << "FOUND SEQEUNCE starting at " << maxCommonValue << endl;
	}
	else
	{
		cout << "... not found\n";
	}



	return maxCommonValue;
}

int Sequence::GetHashCollisionTheoretical(unsigned int numberOfItemsToStore, unsigned int hashTableSize)
{
//	http://en.wikipedia.org/wiki/Birthday_problem#Generalizations
    double intermediateResult =( (double) hashTableSize - 1.0 )/ (double)hashTableSize;
    intermediateResult = pow(intermediateResult, (double) numberOfItemsToStore);
    intermediateResult = intermediateResult*hashTableSize - hashTableSize + numberOfItemsToStore;

    int collisionCount = (int)intermediateResult;
    return collisionCount;
}

//double GetHashCollisionExpected(double numberOfitems, double  hashTableSize)
//{
//	return numberOfitems* numberOfitems / hashTableSize;
//
//	// http://preshing.com/20110504/hash-collision-probabilities
//
//	double NumberToFillTable = hashTableSize * log(hashTableSize);
//	return numberOfitems - NumberToFillTable;
//
//	double exponentNumerator = -numberOfitems* (numberOfitems -1.0);
//	double exponent = 0.5 * exponentNumerator / hashTableSize;
//	double expvalue = exp(exponent);
//	return 1.0 - expvalue;
//}
bool Sequence::createHashTable ()
{
		
	Timer timer;
	timer.Start ();

	allocateHashTable();
	//const int hashLength = 10; 
	//int hashTableLength = 1 << (2*10); //this is binary equivalent to 2^10 
	char* pLib = sequenceData;
	unsigned int * pLibUnsignedInt = (unsigned int *)sequenceData;
	hashCollisions = 0;
	for (int seqLoc = 0; seqLoc < (sequenceLength - NucleotiesInHash); seqLoc++)
	{
		int HashVal =GetHashValue(pLib,seqLoc,NucleotiesInHash);
		int HashValAlternate = GetHashValue(pLibUnsignedInt,seqLoc,NucleotiesInHash);

		//nucValue is the index (hashValue) of the sequence. we store the location of it (i) at this index location
		if (hashTable [HashVal]!=0)//hash Collision! 
		{
			hashCollisions ++;
		}
		hashTable [HashVal] = seqLoc; 
	}

	timer.Stop();
	timeToCreateHash = timer.Elapsed();

	return (true); 


}

bool Sequence::InstrumentedcreateHashTable ()
{
	allocateHashTable();
	//const int hashLength = 10; 
	//int hashTableLength = 1 << (2*10); //this is binary equivalent to 2^10 
	char* pLib = sequenceData;
//	unsigned int * pLibUnsignedInt = (unsigned int *)sequenceData;
	hashCollisions = 0;
	for (int seqLoc = 0; seqLoc < (sequenceLength - NucleotiesInHash); seqLoc++)
	{
		int HashVal =GetHashValue(pLib,seqLoc,NucleotiesInHash);
	//	int HashValAlternate = GetHashValue(pLibUnsignedInt,seqLoc,NucleotiesInHash);

		//nucValue is the index (hashValue) of the sequence. we store the location of it (i) at this index location
		if (hashTable [HashVal]!=0)//hash Collision! 
		{
			hashCollisions ++;
		}
		hashTable [HashVal] = seqLoc; 
	}
	cout << "Collisions = " << hashCollisions << " Expected = " << GetHashCollisionTheoretical( sequenceLength, getHashTableSize())   << endl;
	return (true); 
}

int Sequence::GetActualHashCollsions()
{

	std::vector<unsigned int> hashTableVec(getHashTableSize());
	char* pLib = sequenceData;
	hashCollisions = 0;
	for (int seqLoc = 0; seqLoc < (sequenceLength - NucleotiesInHash); seqLoc++)
	{
		int HashVal =GetHashValue(pLib,seqLoc,NucleotiesInHash);
	
		//nucValue is the index (hashValue) of the sequence. we store the location of it (i) at this index location
		if (hashTableVec [HashVal]!=0)//hash Collision! 
		{
			hashCollisions ++;
		}
		hashTableVec [HashVal] = seqLoc; 
	}

	return hashCollisions;


}

using namespace Concurrency;

bool Sequence::createHashTableWithAMP()
{
	
	Timer timer;
	timer.Start ();
	int HashTableLength = getHashTableSize();

	Concurrency::extent<1> libraryExtent (sequenceLength-12);// need to round number of nucleotides in hash to a multiple of 4 
	Concurrency::array_view<const unsigned int, 1> libraryArrayView (sequenceLength/4,(unsigned int*) sequenceData); // the array view handles copying data before the cPU memory and the GPU memory 

	Concurrency::extent<1> libraryHashExtent (HashTableLength);
	allocateHashTable(); 

	//Compare performance for int and unsigned int and const unsigned int
	Concurrency::array_view<unsigned int, 1> libraryHashArrayView (HashTableLength,(unsigned int*) hashTable); // the array view handles copying data before the cPU memory and the GPU memory 


	//Concurrency::extent<1> probeExtent ((Probe.sequenceLength/4)-3);
	//Concurrency::array_view<const unsigned int, 1> probeArrayView (Probe.sequenceLength/4,(unsigned int*) Probe.sequenceData); // the array view handles copying data before the cPU memory and the GPU memory 

	int NucleotidesPerHash = NucleotiesInHash;
	libraryHashArrayView.discard_data();


	parallel_for_each(libraryExtent, [=] (index<1> idx) restrict (amp) {
	int libraryOffset = idx [0];  // current threads index from libraryExtent, zero to library size 

	///create hash for libraryOffset
	
	  // Calculate the hash for first pattern_length characters of the record (exclude null terminating characters)
// Copyright (c) Microsoft Corp. 
//
// File: ampstrstr.cpp
//
// Licensed under the Apache License, Version 2.0 (the "License"); you may not use this 
// file except in compliance with the License. You may obtain a copy of the License at 
// http://www.apache.org/licenses/LICENSE-2.0  
//  
// THIS CODE IS PROVIDED *AS IS* BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, 
// EITHER EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION ANY IMPLIED WARRANTIES OR 
// CONDITIONS OF TITLE, FITNESS FOR A PARTICULAR PURPOSE, MERCHANTABLITY OR NON-INFRINGEMENT. 
//  
		//This hash function assumes that the characters are encoded four each into an int array
	//  Offest is still in chars
	//  int1   int2  int3
	//[ ACGT] [CGTA] [GTAC]

	//First, we need to figure out how walk bytes in the stream starting at offset
	// Divide by 4 will get the correct integer -- so for example, if we start at byte 2, 2/4 = 0 in integer division
	// we use the modulo (remainder to figure out how much to shift
	

	
		int HashVal =0;
		for (int j = 0; j < NucleotidesPerHash ; j++) //Cannot use Library.NucleotiesInHash because Library object uses Chars
		{
		// Divide by 4 will get the correct integer -- so for example, if we start at byte 2, 2/4 = 0 in integer division
		unsigned int cval =  libraryArrayView ((libraryOffset + j)/4);
		// we use the modulo (remainder to figure out how much to shift
		int remainder = (libraryOffset + j) % 4;

		//Integers are stored in backwards -- least significant byte first  (ACGT in memory corresponds to TGCA as a 32 bit integer)
		// if the remainder is 0, we are at 0, 4, 8, etc    ACGT so we want A
		//			which encodes least significant byte  so we need to shift to right 3*8 bits
		// if the remainder is 1, we are at 1, 5, 9, etc    ACGT so we want C,
		//			so we need to shift to right 1*8 bits
		// if the remainder is 2, we are at 2, 6, 10, etc   ACGT so we want G, so we need to shift to right 2*8 bits
		// if the remainder is 3, we are at 3, 7, 11, etc   ACGT so we want T, so we need to shift to right 3*8 bits
		cval = cval >>  remainder *8;
	
		//now we need to mask off the upper three bytes by ANDing with 0x000000FF;
		cval =cval & 0x000000FF;

	
			int NucValue = 0; 
			switch (cval)
			{
			case 'A': 
				{
					NucValue = 0;
				}
				break;
			case 'G': 
				{
					NucValue = 1;
				}
				break;
			case 'T': 
				{
					NucValue = 2;
				}
				break;
			case 'C': 
				{
					NucValue = 3;
				}
				break;
			default:
				break;
			}

			HashVal = HashVal << 2 ; //The '<<' operator will shift bits to the left when the char on the left is a number and not a 'cout' 
			HashVal += NucValue;
		}//end of Jloop

		

       
		libraryHashArrayView (HashVal) = libraryOffset; 
	///end create hash 


	
	});//end of parallelforeach 
	libraryHashArrayView.synchronize (); 
	timer.Stop();
	timeToCreateHash = timer.Elapsed();
	expectedHashCollisions = GetHashCollisionTheoretical(libraryExtent.size(),getHashTableSize());

	return true;
}

void Sequence::RandomizeSequence(int seed)
{
	PopulateRandomSequence(seed,sequenceLength, sequenceData);
}
	

void Sequence::RetrieveSubseqence(unsigned int start, int length,Sequence & probe)
{
	//create a sequence equal to the desired subsequence length
	probe.init(length,NucleotiesInHash);
	probe.sequenceData = (char*) malloc(length+1);
	memset(probe.sequenceData,0, length+1);
	memcpy(probe.sequenceData,&sequenceData[start],length);
	strcpy(probe.Title , Title);
	strcat(probe.Title, "-subseq@");
	char buffer50[50];
	itoa(start, buffer50,10);
	strcat(probe.Title, buffer50);
	strcat(probe.Title, "-len=");
	itoa(length, buffer50,10);
	strcat(probe.Title, buffer50);
	
}
